Hoefer et al 2023 - Sampling Methods in Reptiles
Sebastian Hoefer
2023-05-27
1 Load packages
Load packages needed for data manipulation and analyses.
library(knitr) #for markdown knitting
library(tidyverse) #for data wrangling etc
library(colorspace) #for working with colours in visualisations
library(cmdstanr) #for cmdstan
library(brms) #for fitting models in STAN
library(rstan) #for interfacing with STAN
library(ggeffects) #for partial effects plots
library(DHARMa) #for residual diagnostics
library(emmeans) #for marginal means etc
library(tidybayes) #for more tidying outputs
library(ggridges) #for creating ridge plots
library(vegan) #for analysing ecological data
library(EcolUtils) #for analysing ecological data
library(patchwork) #for multi-panel figures
library(HDInterval) #for HPD intervals
library(report)2 Functions
Custom ggplot theme for visualisation.
my.theme <- function(){
theme_classic() +
theme(text = element_text(family = "Avenir Next"),
axis.title.y = element_text(margin = margin(t = 0,r = 20,b = 0,l = 0)),
axis.title.x = element_text(margin = margin(t = 20,r = 0,b = 0,l = 0)),
plot.margin = unit(c(5, 10, 5, 10), units = "mm"),
strip.background = element_rect(fill = "#CCCCFF"),
strip.text.x = element_text(size = 20),
axis.title = element_text(size = 20),
axis.text = element_text(size = 18),
legend.text = element_text(size = 15),
legend.title = element_text(size = 15))
}3 Load data
reptile_data <- read.csv("reptile_data.csv")4 Species Richness
4.1 Prepare data
Created an object with coordinates for each plot to assess the effect of latitude.
site <- c(rep("Duval",4), rep("Mourachan",4), rep("Tarcutta",4), rep("Wambiana",4), rep("Undara",4), rep("Rinyirru",4))
site.plot <- c("Duval.DryA","Duval.DryB","Duval.WetA", "Duval.WetB", "Mourachan.DryA", "Mourachan.DryB", "Mourachan.WetA", "Mourachan.WetB", "Tarcutta.DryA", "Tarcutta.DryB", "Tarcutta.WetA", "Tarcutta.WetB", "Wambiana.DryA", "Wambiana.DryB", "Wambiana.WetA", "Wambiana.WetB", "Undara.DryA", "Undara.DryB", "Undara.WetA", "Undara.WetB", "Rinyirru.DryA", "Rinyirru.DryB", "Rinyirru.WetA", "Rinyirru.WetB")
lat <- c(-30.41734901, -30.40221297, -30.41803398, -30.40110799, -27.77898598, -27.77996097, -27.78425401, -27.77876303, -35.36852497, -35.37876002, -35.36012203, -35.36613797, -20.53070999, -20.526226, -20.53458997, -20.53051997, -18.18705803, -18.24565097, -18.18492902, -18.26552999, -15.05447703, -15.04400703, -15.04891698, -15.04438899)
long <- c(151.622667, 151.624706, 151.61276, 151.629483, 149.032006, 148.981, 149.021332, 148.970948, 147.696893, 147.703341, 147.697579, 147.707788, 146.113426, 146.110754, 146.103876, 146.101471, 144.539753, 144.553799, 144.5324, 144.556613, 144.256419, 144.242852, 144.261894, 144.26052)
wet.dry <- rep(c("dry", "dry", "wet", "wet"),6)
coord <- data.frame(site, site.plot, lat, long, wet.dry)Summarised the total number of species per method for each plot across all sites.
reptile_summary <- reptile_data %>%
unite(site.plot, c(site, plot), sep = ".", remove = FALSE) %>%
select(-plot) %>%
group_by(site, site.plot, assessment.method) %>%
summarise(richness = length(unique(scientific.name))) %>%
ungroup() %>%
add_column(wet.dry = ifelse(.$site.plot %in% c("Duval.WetA", "Duval.WetB",
"Tarcutta.WetA", "Tarcutta.WetB",
"Mourachan.WetA", "Mourachan.WetB",
"Wambiana.WetA", "Wambiana.WetB",
"Undara.WetA", "Undara.WetB",
"Rinyirru.WetA", "Rinyirru.WetB"), "wet", "dry"))The data frame was missing zero values for methods that didn’t capture any species at a given plot. The following adds zeros for those instances.
# Created object with all combinations of plots and methods to include zero values
plots <- unique(reptile_summary$site.plot)
methods <- unique(reptile_summary$assessment.method)
zeros <- crossing(plots, methods)
# Added reference column for future join
zeros$plots.methods <- paste(zeros$plots, zeros$methods)
# Selected only columns relevant to join
reptile_summary2 <- reptile_summary %>%
select("site.plot", "assessment.method", "richness")
# Added reference column for future join
reptile_summary2$plots.methods <- paste(reptile_summary2$site.plot, reptile_summary2$assessment.method)
# Joined objects to include zero values for methods across plots
reptile_summary_all <- merge(reptile_summary2, zeros, by = "plots.methods",
all.x = T, all.y = T) %>%
select("plots", "methods", "richness") %>%
replace(is.na(.), 0) %>%
rename(site.plot = plots, assessment.method = methods) %>%
left_join(coord, by = "site.plot") %>%
mutate(site = factor(site, levels = c("Tarcutta", "Duval", "Mourachan", "Wambiana",
"Undara", "Rinyirru")),
site.plot = factor(site.plot, levels = c("Tarcutta.DryA", "Tarcutta.DryB", "Tarcutta.WetA", "Tarcutta.WetB", "Duval.DryA", "Duval.DryB", "Duval.WetA", "Duval.WetB", "Mourachan.DryA", "Mourachan.DryB", "Mourachan.WetA", "Mourachan.WetB","Wambiana.DryA", "Wambiana.DryB", "Wambiana.WetA", "Wambiana.WetB", "Undara.DryA", "Undara.DryB", "Undara.WetA", "Undara.WetB", "Rinyirru.DryA", "Rinyirru.DryB", "Rinyirru.WetA", "Rinyirru.WetB")),
assessment.method = factor(assessment.method, levels = c("pitfall", "funnel","incidentals","spotlighting","cover board","camera")),
wet.dry = factor(wet.dry),
richness = as.numeric(richness))4.2 Bayesian Analysis
4.2.1 Priors
Defining priors for Bayesian models.
Priors for the intercept.
reptile_summary_all %>% summarise(median(log(richness)))## median(log(richness))
## 1 1.098612#1.1
reptile_summary_all %>% summarise(mad(log(richness)))## mad(log(richness))
## 1 1.454177#1.5Priors for the slope.
log(sd(reptile_summary_all$richness))/apply(model.matrix(~assessment.method*lat, data = reptile_summary_all), 2, sd)## (Intercept) assessment.methodfunnel
## Inf 3.6955179
## assessment.methodincidentals assessment.methodspotlighting
## 3.6955179 3.6955179
## assessment.methodcover board assessment.methodcamera
## 3.6955179 3.6955179
## lat assessment.methodfunnel:lat
## 0.1921267 0.1433234
## assessment.methodincidentals:lat assessment.methodspotlighting:lat
## 0.1433234 0.1433234
## assessment.methodcover board:lat assessment.methodcamera:lat
## 0.1433234 0.1433234#3.74.2.2 Fit model
4.2.2.1 Model 1
priors1 <- prior(normal(1.1,1.5), class = "Intercept") +
prior(normal(0,3.7), class = "b") +
prior(student_t(3,0,3.7), class = "sd")
methods.form1 <- bf(richness ~ assessment.method*scale(lat, scale = FALSE) + (1|site/site.plot),
family = poisson(link = "log"))
methods.brm1 <- brm(methods.form1,
data = reptile_summary_all,
prior = priors1,
sample_prior = "only",
refresh = 0,
chains = 3, cores = 3,
iter = 5000,
thin = 5,
seed = 8,
warmup = 1000,
backend = 'cmdstanr',
save_pars = save_pars(all = TRUE))## Running MCMC with 3 parallel chains...
##
## Chain 1 finished in 0.2 seconds.
## Chain 2 finished in 0.2 seconds.
## Chain 3 finished in 0.2 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.3 seconds.4.2.2.1.1 Predictions
methods.brm1 %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)
methods.brm1 %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)
4.2.2.1.2 Prior vs Posterior
methods.brm1b <- methods.brm1 %>%
update(sample_prior = "yes",
refresh = 0,
seed = 8,
cores = 3,
backend = "cmdstanr",
save_pars = save_pars(all = TRUE))## Running MCMC with 3 parallel chains...
##
## Chain 1 finished in 2.1 seconds.
## Chain 3 finished in 2.2 seconds.
## Chain 2 finished in 3.3 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 2.5 seconds.
## Total execution time: 3.4 seconds.4.2.2.1.3 Predictions
methods.brm1b %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)
methods.brm1b %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)
4.2.2.2 Model 2
priors2 <- prior(normal(1.1,1.5), class = "Intercept") +
prior(normal(0,3.7), class = "b") +
prior(student_t(3,0,3.7), class = "sd")
methods.form2 <- bf(richness ~ assessment.method*scale(lat, scale = FALSE) + (1|site/site.plot),
family="negbinomial")
methods.brm2 <- brm(methods.form2,
data = reptile_summary_all,
prior = priors2,
sample_prior = "only",
refresh = 0,
chains = 3, cores = 3,
iter = 5000,
thin = 5,
seed = 8,
warmup = 1000,
backend = 'cmdstanr',
save_pars = save_pars(all = TRUE))## Running MCMC with 3 parallel chains...
##
## Chain 2 finished in 1.9 seconds.
## Chain 3 finished in 3.5 seconds.
## Chain 1 finished in 3.9 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 3.1 seconds.
## Total execution time: 4.0 seconds.4.2.2.2.1 Predictions
methods.brm2 %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)
methods.brm2 %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)
4.2.2.2.2 Prior vs Posterior
methods.brm2b <- methods.brm2 %>%
update(sample_prior = "yes",
refresh = 0,
seed = 8,
cores = 3,
backend = "cmdstanr",
save_pars = save_pars(all = TRUE))## Running MCMC with 3 parallel chains...
##
## Chain 1 finished in 2.9 seconds.
## Chain 3 finished in 2.9 seconds.
## Chain 2 finished in 3.0 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 3.0 seconds.
## Total execution time: 3.1 seconds.4.2.2.2.3 Predictions
methods.brm2b %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)
methods.brm2b %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)
4.2.2.3 Model 3
priors3 <- prior(normal(1.1,1.5), class = "Intercept") +
prior(normal(0,3.7), class = "b") +
prior(student_t(3,0,3.7), class = "sd")
methods.form3 <- bf(richness ~ assessment.method*scale(lat, scale = FALSE) + (1|site/site.plot),
family="negbinomial2")
methods.brm3 <- brm(methods.form3,
data = reptile_summary_all,
prior = priors3,
sample_prior = "only",
refresh = 0,
chains = 3, cores = 3,
iter = 5000,
thin = 5,
seed = 1,
warmup = 1000,
control = list(adapt_delta=0.99),
backend = 'cmdstanr',
save_pars = save_pars(all = TRUE))## Running MCMC with 3 parallel chains...
##
## Chain 2 finished in 0.5 seconds.
## Chain 1 finished in 0.6 seconds.
## Chain 3 finished in 0.5 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 0.5 seconds.
## Total execution time: 0.7 seconds.4.2.2.3.1 Predictions
methods.brm3 %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)
methods.brm3 %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)
4.2.2.3.2 Prior vs Posterior
methods.brm3b <- methods.brm3 %>%
update(sample_prior = "yes",
refresh = 0,
seed = 2,
cores = 3,
control = list(adapt_delta=0.99),
backend = "cmdstanr",
save_pars = save_pars(all = TRUE))## Running MCMC with 3 parallel chains...
##
## Chain 1 finished in 6.6 seconds.
## Chain 2 finished in 6.6 seconds.
## Chain 3 finished in 9.5 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 7.5 seconds.
## Total execution time: 9.5 seconds.4.2.2.3.3 Predictions
methods.brm3b %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)
methods.brm3b %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)
4.2.2.4 Compare Models
(l.1b <- methods.brm1b %>% loo())##
## Computed from 2400 by 144 log-likelihood matrix
##
## Estimate SE
## elpd_loo -265.3 11.2
## p_loo 16.5 2.2
## looic 530.6 22.3
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 142 98.6% 427
## (0.5, 0.7] (ok) 2 1.4% 759
## (0.7, 1] (bad) 0 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
##
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.(l.2b <- methods.brm2b %>% loo())##
## Computed from 2400 by 144 log-likelihood matrix
##
## Estimate SE
## elpd_loo -266.7 11.1
## p_loo 16.0 2.1
## looic 533.5 22.2
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 143 99.3% 456
## (0.5, 0.7] (ok) 1 0.7% 660
## (0.7, 1] (bad) 0 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
##
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.(l.3b <- methods.brm3b %>% loo())##
## Computed from 2400 by 144 log-likelihood matrix
##
## Estimate SE
## elpd_loo -266.8 11.2
## p_loo 16.2 2.2
## looic 533.5 22.3
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 141 97.9% 642
## (0.5, 0.7] (ok) 3 2.1% 164
## (0.7, 1] (bad) 0 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
##
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.loo_compare(loo(methods.brm1b), loo(methods.brm2b), loo(methods.brm3b))## elpd_diff se_diff
## methods.brm1b 0.0 0.0
## methods.brm2b -1.4 0.4
## methods.brm3b -1.5 0.4Model methods.brm1b was selected as best model based on
loo estimates.
4.2.3 Diagnostics
methods.brm1b %>% hypothesis("scalelatscaleEQFALSE = 0") %>% plot
4.2.3.1 Trace plots
methods.brm1b$fit %>% stan_trace()
#### Autocorrelation plots
methods.brm1b$fit %>% stan_ac()
#### Rhat statistic
methods.brm1b$fit %>% stan_rhat()
#### Effective sampling size
methods.brm1b$fit %>% stan_ess()
#### Posterior predictive check plot
methods.brm1b %>% pp_check(type = "dens_overlay", ndraws = 200)
#### DHARMa residuals
set.seed(6)
preds <- posterior_predict(methods.brm1b, ndraws=250, summary=FALSE)
method.resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = reptile_summary_all$richness,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE)
method.resids %>% plot()
#### Dispersion test
method.resids %>% testDispersion()
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.57235, p-value < 2.2e-16
## alternative hypothesis: two.sided4.2.3.2 Zero-inflation test
method.resids %>% testZeroInflation()
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.1598, p-value = 0.4
## alternative hypothesis: two.sided4.2.4 Investigation
4.2.4.1 Methods pairwise comparison for each site
newdata_lat <- with(reptile_summary_all,list(lat = c(-35.37876,-30.40221,-27.77876,
-20.53459,-18.26553,-15.04439)))
(diff.methods3 <- methods.brm1b %>%
emmeans(~assessment.method|lat, at = newdata_lat) %>%
pairs() %>%
gather_emmeans_draws() %>%
mutate(f.change = exp(.value)) %>%
summarise("Average fractional change" = median(f.change),
"Lower HDI" = HDInterval::hdi(f.change)[1],
"Upper HDI" = HDInterval::hdi(f.change)[2],
"Probability of difference" = sum(.value > 0)/n()) %>% #to see if there is any change
arrange(lat) %>%
rename("Site" = lat) %>%
mutate(Site = as.factor(Site)) %>%
mutate(Site = fct_recode(Site, "Tarcutta" = '-35.37876', "Duval" = '-30.40221',
"Mourachan" = '-27.77876', "Wambiana" = '-20.53459',
"Undara" = '-18.26553', "Rinyirru" = '-15.04439')))## # A tibble: 90 × 6
## # Groups: contrast [15]
## contrast Site Average fractional c…¹ `Lower HDI` `Upper HDI`
## <fct> <fct> <dbl> <dbl> <dbl>
## 1 pitfall - funnel Tarc… 1.05 0.615 1.71
## 2 pitfall - incidentals Tarc… 4.04 1.64 7.72
## 3 pitfall - spotlighting Tarc… 4.81 2.03 9.65
## 4 pitfall - cover board Tarc… 8.03 2.53 18.7
## 5 pitfall - camera Tarc… 8.42 1.79 24.0
## 6 funnel - incidentals Tarc… 3.79 1.59 7.34
## 7 funnel - spotlighting Tarc… 4.55 1.67 8.83
## 8 funnel - cover board Tarc… 7.58 2.35 17.8
## 9 funnel - camera Tarc… 8.03 1.96 23.6
## 10 incidentals - spotlight… Tarc… 1.19 0.356 2.65
## # ℹ 80 more rows
## # ℹ abbreviated name: ¹`Average fractional change`
## # ℹ 1 more variable: `Probability of difference` <dbl>4.2.4.2 Methods pairwise comparison for the average latitude
(diff.meth.avg <- methods.brm1b %>%
emmeans(~assessment.method|lat) %>%
pairs() %>%
gather_emmeans_draws() %>%
mutate(f.change = exp(.value)) %>%
summarise("Average fractional change" = median(f.change),
"Lower HDI" = HDInterval::hdi(f.change)[1],
"Upper HDI" = HDInterval::hdi(f.change)[2],
"Probability of difference" = sum(.value > 0)/n()) %>%
select(-lat))## # A tibble: 15 × 5
## # Groups: contrast [15]
## contrast Average fractional chang…¹ `Lower HDI` `Upper HDI`
## <fct> <dbl> <dbl> <dbl>
## 1 pitfall - funnel 1.02 0.801 1.28
## 2 pitfall - incidentals 2.50 1.77 3.37
## 3 pitfall - spotlighting 3.03 2.12 4.21
## 4 pitfall - cover board 4.80 2.98 7.09
## 5 pitfall - camera 11.8 6.74 19.8
## 6 funnel - incidentals 2.45 1.79 3.40
## 7 funnel - spotlighting 2.96 2.04 4.11
## 8 funnel - cover board 4.67 2.78 6.79
## 9 funnel - camera 11.5 5.89 19.0
## 10 incidentals - spotlighting 1.21 0.776 1.75
## 11 incidentals - cover board 1.91 1.13 2.98
## 12 incidentals - camera 4.69 2.25 7.97
## 13 spotlighting - cover board 1.58 0.945 2.51
## 14 spotlighting - camera 3.87 1.87 6.66
## 15 cover board - camera 2.47 1.08 4.40
## # ℹ abbreviated name: ¹`Average fractional change`
## # ℹ 1 more variable: `Probability of difference` <dbl>4.2.5 Visualisation
Fractional change
methods.em <- methods.brm1b %>%
emmeans(~assessment.method|lat) %>%
pairs() %>%
gather_emmeans_draws() %>%
mutate(rate = exp(.value))
(methods.mag <- methods.em %>%
ggplot() +
geom_density_ridges_gradient(aes(x=rate, y=fct_reorder(contrast,rate)),
alpha = 0.4, col = "white",
quantile_lines = TRUE, quantiles = c(0.025, 0.975),
show.legend = FALSE, fill = "#05596E") +
geom_vline(xintercept = 1, linetype = "dashed") +
scale_x_continuous("Fractional change", trans = "log2",
breaks = c(1, 2, 5, 10, 50)) +
scale_y_discrete(name = "") +
my.theme())
Proportion of richness across latitudes
lat.list <- with(reptile_summary_all, list(lat = seq(min(lat), max(lat), length = 100)))
newdata <- emmeans(methods.brm1b, ~lat|assessment.method, at=lat.list, type = "response") %>%
gather_emmeans_draws() %>%
mutate(richness = exp(.value)) %>%
group_by(lat, .draw) %>%
mutate(sum_rate = sum(richness)) %>%
ungroup() %>%
mutate(proportion = richness/sum_rate) %>%
filter(.draw %in% sample(1:2400, 200))
newdata3 <- emmeans(methods.brm1b, ~lat|assessment.method, at=lat.list, type = "response") %>%
as.data.frame %>%
group_by(lat) %>%
mutate(sum_rate.avg = sum(rate),
proportion.avg = rate/sum_rate.avg) %>%
ungroup()
reptile_summary_all_prop <- reptile_summary_all %>%
group_by(lat) %>%
mutate(sum_rich = sum(richness),
proportion = richness/sum_rich)
(methods.spag <- newdata %>%
ggplot() +
geom_jitter(data = reptile_summary_all_prop, aes(x = lat, y = proportion, col = assessment.method),
height = 0, width = 0.2, alpha = 0.4) +
geom_line(aes(lat, proportion, col = assessment.method, group = interaction(assessment.method,.draw)), alpha=0.07) +
geom_line(data = newdata3, aes(lat, proportion.avg, group = assessment.method), linewidth = 1.7, col = "black") +
geom_line(data = newdata3, aes(lat, proportion.avg, col = assessment.method), linewidth = 1) +
my.theme() +
scale_y_continuous(name = "Proportion of total species richness",
labels = scales::percent,
breaks = seq(0, 0.5, by = 0.1),
limits = c(-0.025,0.5)) +
scale_x_continuous(name = "", breaks=c(-35,-30,-25,-20,-15),
labels=c("35ºS", "30ºS", "25ºS", "20ºS", "15ºS")) +
scale_fill_manual(values = c("#FF8300", "#D14103","#0CC170","black","#4E84C4","#8348D8"),
name = "Sampling Methods", labels = c("Pitfall Trap","Funnel Trap","Incidental","Spotlighting","Arboreal Cover Board","Camera Trap")) +
scale_colour_manual(values = c("#FF8300", "#D14103","#0CC170","black","#4E84C4","#8348D8"),
name = "", labels = c("Pitfall Trap","Funnel Trap","Incidental","Spotlighting","Arboreal Cover Board","Camera Trap")) +
theme(legend.text = element_text(size = 16, margin = margin(t = 5)),
legend.position = c(y=0.45, x=-0.07),
legend.background = element_rect(fill = "transparent")) +
guides(colour = guide_legend(nrow = 1)) +
annotate("text", x = c(-35.37876,-30.40221,-27.77876,-20.53459,-18.26553,-15.04439),
y = c(rep(-0.025, 6)),
label = c("Tarcutta","Duval","Mourachan","Wambiana","Undara","Rinyirru"), size = 5))
5 Survey Effort
Species accumulation curves
5.1 Prepare data
Merged duplicate rows for duplicate observations and transformed the dataset to a wide format.
reptile_community <- reptile_data %>%
select(site, date, assessment.method,scientific.name) %>% #selected relevant grouping variables
group_by(across(everything())) %>% #groups by all available columns
mutate(number = n()) %>% #created a numbers column that counts the number of duplicate observations
ungroup() %>%
unique() %>% #merged the duplicate rows
pivot_wider(names_from = "scientific.name", #transformed dataset from long to wide format
values_from = "number",
values_fill = list(number=0)) %>% #replaced NAs with 0
arrange(site)Add a category for the total richness.
result.i <- vector("list",length(unique(reptile_community$site))) #created an empty list that was filled as the loop ran
for(i in 1:length(unique(reptile_community$site))){ #looped over each study site
site.i <- reptile_community[reptile_community$site == unique(reptile_community$site)[i],] #subset to the site
result.a <- vector("list",length(unique(site.i$date))) #created an empty list to be filled by each date
for(a in 1:length(unique(site.i$date))){ #looped over each date
date.a <- site.i[site.i$date == unique(site.i$date)[a],] #subset to the date
meta.a <- cbind.data.frame(date.a$site[1],date.a$date[1],c("total"))
date.a <- date.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
total.a <- cbind.data.frame(meta.a,rbind.data.frame(colSums(date.a))) #calculated column sums
colnames(total.a) <- colnames(reptile_community)
result.a[[a]] <- total.a} #added method results to list
result.i[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop
result.i <- do.call("rbind.data.frame",result.i)
reptile_community2 <- rbind.data.frame(reptile_community, result.i)Add a category for the combined pitfall and funnel richness.
result.j <- vector("list",length(unique(reptile_community$site))) #created an empty list that was filled as the loop ran
for(i in 1:length(unique(reptile_community$site))){ #looped over each study site
site.i <- reptile_community[reptile_community$site == unique(reptile_community$site)[i],] #subset to the site
result.a <- vector("list",length(unique(site.i$date))) #created an empty list to be filled by each date
for(a in 1:length(unique(site.i$date))){ #looped over each date
date.a <- site.i[site.i$date == unique(site.i$date)[a],] %>%
filter(assessment.method %in% c("funnel", "pitfall"))
meta.a <- cbind.data.frame(date.a$site[1],date.a$date[1],c("pitfunnel"))
date.a <- date.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
total.a <- cbind.data.frame(meta.a,rbind.data.frame(colSums(date.a))) #calculated column sums
colnames(total.a) <- colnames(reptile_community)
result.a[[a]] <- total.a} #added method results to list
result.j[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop
result.j <- do.call("rbind.data.frame",result.j)
reptile_community2 <- rbind.data.frame(reptile_community2, result.j) %>%
drop_na()5.2 Loop
Created a loop that applied specaccum() to all methods across all sites. The output were values for richness of each method at each site over the course of the survey.
# Created a data frame of zeros with 30 rows and the same number of columns as the total number of species.
extra.rows <- data.frame(matrix(NA, nrow = 30, ncol = (ncol(reptile_community2)))) %>%
replace(is.na(.), 0) %>%
set_names(names(reptile_community2)) %>%
select(-(1:3))
result.i <- vector("list",length(unique(reptile_community2$site))) #created an empty list that was filled as the loop ran
set.seed(8) #set seed for reproducibility
for(i in 1:length(unique(reptile_community2$site))){ #looped over each study site
site.i <- reptile_community2[reptile_community2$site == unique(reptile_community2$site)[i],] #subset to the site
dates.i <- ifelse(unique(site.i$site) %in% c("Tarcutta", "Undara", "Wambiana", "Rinyirru"), 28, 21) #manual addition of dates depending on site to avoid any issues from 0 animals being detected by any method on a given day
result.a <- vector("list",length(unique(site.i$assessment.method))) #created an empty list to be filled by each method
for(a in 1:length(unique(site.i$assessment.method))){ #looped over each method
method.a <- site.i[site.i$assessment.method == unique(site.i$assessment.method)[a],] #subset to the method
if(nrow(method.a) > 1){ #only proceeded if detections occurred on more than 1 day for specaccum() to work properly
sampling.a <- method.a$assessment.method[1] #saved the name of the method
method.a <- method.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
if(nrow(method.a)!=dates.i){ #if there were missing days (i.e., days when nothing was found, use the extra.rows object to add rows of zeros)
method.a <- rbind.data.frame(method.a,extra.rows[1:(dates.i-nrow(method.a)),])}
accum.a <- specaccum(method.a,"random") #calculated data
#created data frame of metadata, richness, and standard deviation
accum.a <- cbind.data.frame(rep(site.i$site[1],nrow(method.a)),rep(sampling.a,nrow(method.a)),accum.a$richness,accum.a$sd, 1:nrow(method.a))
colnames(accum.a) <- c("site","assessment.method","richness","sd", "day")
accum.a$lower.ci <- accum.a$richness-qnorm(0.975)*accum.a$sd/sqrt(100) #calculated lower 95% CI
accum.a$upper.ci <- accum.a$richness+qnorm(0.975)*accum.a$sd/sqrt(100) #calculated upper 95% CI
result.a[[a]] <- accum.a}} #added method results to list
result.i[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loopTransformed resulting list to a data frame.
reptiles_wide.result <- do.call("rbind.data.frame",result.i) %>% #compressed result.i list into data frame
mutate(assessment.method = factor(assessment.method, levels = c("pitfall",
"funnel",
"incidentals",
"spotlighting",
"cover board",
"camera",
"total",
"pitfunnel")),
site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
"Duval", "Tarcutta")))5.3 Visualisation
Plot of rarefaction curves for each method split by sites.
(specaccum_reptiles <- ggplot(reptiles_wide.result, aes(x = day, y = richness, col = assessment.method)) +
geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci,fill=after_scale(alpha(colour, 0.3))),
linetype = 0.1) +
geom_vline(xintercept=c(7,14,21, 28), linetype="dotted", colour = "black") +
geom_line(aes(group = assessment.method)) +
my.theme() +
facet_wrap(~site) +
theme(panel.border = element_rect(fill = NA, color = "black", linewidth = 1.5),
strip.background = element_rect(fill = "lightgrey"),
axis.title = element_text(size = 26),
axis.text = element_text(size = 26),
legend.text = element_text(size = 26),
legend.position = "bottom",
strip.text.x = element_text(size = 26)) +
scale_y_continuous(expand = expansion(mult = c(0, 0.07)), breaks = seq(0, 45, by = 5),
name = "Species Richness") +
scale_x_continuous(breaks = seq(0, 28, by = 7),
name = "Survey Effort (Days)") +
scale_colour_manual(values = c("#FF8300","#D14103","#0CC170","black",
"#4E84C4","#8348D8","#FBA6EA","#E8BB5A"),
name = "",
labels = c("Pitfall Trap","Funnel Trap","Incidental",
"Spotlighting","Arboreal Cover Board", "Camera Trap",
"Total Richness", "Pit + Funnel")))
6 Community composition
6.1 Prepare data
reptile_community3 <- reptile_data %>%
filter(!recapture %in% c("y")) %>%
select(site, plot, assessment.method,scientific.name) %>% #selected relevant grouping variables
group_by(across(everything())) %>% #groups by all available columns
mutate(number = n()) %>% #created a numbers column that counts the number of duplicate observations
ungroup() %>%
unique() %>% #merged the duplicate rows
pivot_wider(names_from = "scientific.name", #transformed dataset from long to wide format
values_from = "number",
values_fill = list(number=0)) %>% #replaced NAs with 0
arrange(site) %>%
mutate(assessment.method = factor(assessment.method))duv <- reptile_community3 %>% filter(site=="Duval") %>% select(where(~ any(. != 0)))
mou <- reptile_community3 %>% filter(site=="Mourachan") %>% select(where(~ any(. != 0)))
rin <- reptile_community3 %>% filter(site=="Rinyirru") %>% select(where(~ any(. != 0)))
tar <- reptile_community3 %>% filter(site=="Tarcutta") %>% select(where(~ any(. != 0)))
und <- reptile_community3 %>% filter(site=="Undara") %>% select(where(~ any(. != 0)))
wam <- reptile_community3 %>% filter(site=="Wambiana") %>% select(where(~ any(. != 0)))
all <- reptile_community3 %>% select(where(~ any(. != 0)), -"Chelodina longicollis") %>% filter(rowSums(.[,4:ncol(.)]) > 0)#Function decostand() standardises proportions (method = "total") by rows (MARGIN = 1)
rep_com_bray_duv <- duv %>%
select(where(is.numeric)) %>%
decostand(method="total", MARGIN = 1)
rep_com_bray_mou <- mou %>%
select(where(is.numeric)) %>%
decostand(method="total", MARGIN = 1)
rep_com_bray_rin <- rin %>%
select(where(is.numeric)) %>%
decostand(method="total", MARGIN = 1)
rep_com_bray_tar <- tar %>%
select(where(is.numeric)) %>%
decostand(method="total", MARGIN = 1)
rep_com_bray_und <- und %>%
select(where(is.numeric)) %>%
decostand(method="total", MARGIN = 1)
rep_com_bray_wam <- wam %>%
select(where(is.numeric)) %>%
decostand(method="total", MARGIN = 1)
rep_com_bray_all <- all %>%
select(where(is.numeric)) %>%
decostand(method="total", MARGIN = 1)#Function decostand() standardises presence/absence (method = "pa") by rows (MARGIN = 1)
rep_com_jacc_duv <- duv %>%
select(where(is.numeric)) %>%
decostand(method="pa", MARGIN = 1)
rep_com_jacc_mou <- mou %>%
select(where(is.numeric)) %>%
decostand(method="pa", MARGIN = 1)
rep_com_jacc_rin <- rin %>%
select(where(is.numeric)) %>%
decostand(method="pa", MARGIN = 1)
rep_com_jacc_tar <- tar %>%
select(where(is.numeric)) %>%
decostand(method="pa", MARGIN = 1)
rep_com_jacc_und <- und %>%
select(where(is.numeric)) %>%
decostand(method="pa", MARGIN = 1)
rep_com_jacc_wam <- wam %>%
select(where(is.numeric)) %>%
decostand(method="pa", MARGIN = 1)
rep_com_jacc_all <- all %>%
select(where(is.numeric)) %>%
decostand(method="pa", MARGIN = 1)6.2 NMDS Ordination
6.2.1 Bray-Curtis dissimilarity
set.seed(1)
all_nmds_bray <- metaMDS(rep_com_bray_all, distance="bray", k=3,trymax=100, noshare=0.1)
rin_nmds_bray <- metaMDS(rep_com_bray_rin, distance="bray", k=2,trymax=100)
und_nmds_bray <- metaMDS(rep_com_bray_und, distance="bray", k=2,trymax=100)
wam_nmds_bray <- metaMDS(rep_com_bray_wam, distance="bray", k=2,trymax=100)
mou_nmds_bray <- metaMDS(rep_com_bray_mou, distance="bray", k=2,trymax=1000)
#removed due to insufficient data
duv_nmds_bray <- metaMDS(rep_com_bray_duv, distance="bray", k=2,trymax=100)## Warning in metaMDS(rep_com_bray_duv, distance = "bray", k = 2, trymax = 100):
## stress is (nearly) zero: you may have insufficient datatar_nmds_bray <- metaMDS(rep_com_bray_tar, distance="bray", k=2,trymax=100)## Warning in metaMDS(rep_com_bray_tar, distance = "bray", k = 2, trymax = 100):
## stress is (nearly) zero: you may have insufficient data6.2.2 Jaccard dissimilarity
set.seed(11)
all_nmds_jacc <- metaMDS(rep_com_jacc_all, distance="jaccard", k=2,trymax=100, noshare=0.1)
rin_nmds_jacc <- metaMDS(rep_com_jacc_rin, distance="jaccard", k=2,trymax=100)
und_nmds_jacc <- metaMDS(rep_com_jacc_und, distance="jaccard", k=2,trymax=100)
wam_nmds_jacc <- metaMDS(rep_com_jacc_wam, distance="jaccard", k=2,trymax=100)
mou_nmds_jacc <- metaMDS(rep_com_jacc_mou, distance="jaccard", k=2,trymax=100)
#removed due to insufficient data
duv_nmds_jacc <- metaMDS(rep_com_jacc_duv, distance="jaccard", k=2,trymax=100)## Warning in metaMDS(rep_com_jacc_duv, distance = "jaccard", k = 2, trymax =
## 100): stress is (nearly) zero: you may have insufficient datatar_nmds_jacc <- metaMDS(rep_com_jacc_tar, distance="jaccard", k=2,trymax=100)## Warning in metaMDS(rep_com_jacc_tar, distance = "jaccard", k = 2, trymax =
## 100): stress is (nearly) zero: you may have insufficient data6.3 Visualisation
6.3.1 Bray-Curtis dissimilarity
create_plot <- function(data, nmds_data) {
methods <- nmds_data$points
species <- as.data.frame(nmds_data$species) %>% mutate(species = row.names(.))
df <- cbind.data.frame(data[,1:3], methods)
gg <- merge(df, aggregate(cbind(mean.x=MDS1,mean.y=MDS2)~assessment.method, df, mean), by="assessment.method") %>%
mutate(assessment.method = factor(assessment.method, levels = c("pitfall","funnel","incidentals","spotlighting","cover board","camera","total")))
ggplot(gg, aes(MDS1,MDS2,color=assessment.method)) +
geom_point(size=3) +
geom_point(aes(x=mean.x,y=mean.y),size=5) +
scale_colour_manual(values = c("#FF8300", "#D14103","#0CC170","black","#4E84C4","#8348D8"),
name = "",
labels = c("Pitfall Trap","Funnel Trap","Incidental","Spotlighting","Arboreal Cover Board","Camera Trap")) +
geom_segment(aes(x=mean.x, y=mean.y, xend=MDS1, yend=MDS2), alpha = 0.2, linewidth = 1.5) +
my.theme() + theme(legend.position = "bottom")
}
mds.all.plot.bray.centroid <- create_plot(all, all_nmds_bray)
mds.rin.plot.bray.centroid <- create_plot(rin, rin_nmds_bray)
mds.und.plot.bray.centroid <- create_plot(und, und_nmds_bray)
mds.wam.plot.bray.centroid <- create_plot(wam, wam_nmds_bray)
mds.mou.plot.bray.centroid <- create_plot(mou, mou_nmds_bray)
(mds.bray.all <- (mds.all.plot.bray.centroid + ggtitle('All sites') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
(mds.rin.plot.bray.centroid + ggtitle('Rinyirru') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
(mds.und.plot.bray.centroid + ggtitle('Undara') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
(mds.wam.plot.bray.centroid + ggtitle('Wambiana') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
(mds.mou.plot.bray.centroid + ggtitle('Mourachan') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
theme(legend.position = c(2,0.65), legend.text = element_text(size = 24)))
6.3.2 Jaccard dissimilarity
create_plot <- function(data, nmds_data) {
methods <- nmds_data$points
species <- as.data.frame(nmds_data$species) %>% mutate(species = row.names(.))
df <- cbind.data.frame(data[,1:3], methods)
gg <- merge(df, aggregate(cbind(mean.x=MDS1,mean.y=MDS2)~assessment.method, df, mean), by="assessment.method") %>%
mutate(assessment.method = factor(assessment.method, levels = c("pitfall","funnel","incidentals","spotlighting","cover board","camera","total")))
ggplot(gg, aes(MDS1,MDS2,color=assessment.method)) +
geom_point(size=3) +
geom_point(aes(x=mean.x,y=mean.y),size=5) +
scale_colour_manual(values = c("#FF8300", "#D14103","#0CC170","black","#4E84C4","#8348D8"),
name = "",
labels = c("Pitfall Trap","Funnel Trap","Incidental","Spotlighting","Arboreal Cover Board","Camera Trap")) +
geom_segment(aes(x=mean.x, y=mean.y, xend=MDS1, yend=MDS2), alpha = 0.2, linewidth = 1.5) +
my.theme() + theme(legend.position = "bottom")
}
mds.all.plot.jacc.centroid <- create_plot(all, all_nmds_jacc)
mds.rin.plot.jacc.centroid <- create_plot(rin, rin_nmds_jacc)
mds.und.plot.jacc.centroid <- create_plot(und, und_nmds_jacc)
mds.wam.plot.jacc.centroid <- create_plot(wam, wam_nmds_jacc)
mds.mou.plot.jacc.centroid <- create_plot(mou, mou_nmds_jacc)
(mds.jacc.all <- (mds.all.plot.jacc.centroid + ggtitle('All sites') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
(mds.rin.plot.jacc.centroid + ggtitle('Rinyirru') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
(mds.und.plot.jacc.centroid + ggtitle('Undara') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
(mds.wam.plot.jacc.centroid + ggtitle('Wambiana') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
(mds.mou.plot.jacc.centroid + ggtitle('Mourachan') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
theme(legend.position = c(2,0.65), legend.text = element_text(size = 24)))
6.4 Investigation
6.5 Adonis
6.5.1 Bray-Curtis dissimilarity
all_b <- all %>% mutate(assessment.method = factor(assessment.method))
all.dist <- vegdist(all[4:ncol(all)], 'bray')
(all.adonis <- adonis2(all.dist ~ assessment.method, data=all_b))## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = all.dist ~ assessment.method, data = all_b)
## Df SumOfSqs R2 F Pr(>F)
## assessment.method 5 10.309 0.21227 5.5511 0.001 ***
## Residual 103 38.257 0.78773
## Total 108 48.566 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Permutation test: significant difference between assessment methods for reptile communities.
#Let's which method is different from which
(adonis.methods.bray <- adonis.pair(all.dist, all$assessment.method))## combination SumsOfSqs MeanSqs F.Model R2
## 1 camera <-> cover board 2.8879066 2.8879066 9.3629077 0.27247135
## 2 camera <-> funnel 2.7328065 2.7328065 7.5459644 0.19081503
## 3 camera <-> incidentals 2.4035118 2.4035118 6.8296694 0.20803345
## 4 camera <-> pitfall 2.7403269 2.7403269 7.4979872 0.18983213
## 5 camera <-> spotlighting 3.1100425 3.1100425 11.2489845 0.31912932
## 6 cover board <-> funnel 2.6722775 2.6722775 6.9302710 0.15088679
## 7 cover board <-> incidentals 1.5522584 1.5522584 4.0656878 0.10968872
## 8 cover board <-> pitfall 2.2871233 2.2871233 5.8898007 0.13120577
## 9 cover board <-> spotlighting 1.0719637 1.0719637 3.2952071 0.09608360
## 10 funnel <-> incidentals 1.5296767 1.5296767 3.7134502 0.08494983
## 11 funnel <-> pitfall 0.3802544 0.3802544 0.9195332 0.01959809
## 12 funnel <-> spotlighting 2.8043239 2.8043239 7.6321637 0.16725404
## 13 incidentals <-> pitfall 1.4649808 1.4649808 3.5336081 0.08116966
## 14 incidentals <-> spotlighting 1.5851769 1.5851769 4.4019111 0.12092527
## 15 pitfall <-> spotlighting 2.7086276 2.7086276 7.3160470 0.16144495
## P.value P.value.corrected
## 1 0.000999001 0.001152693
## 2 0.000999001 0.001152693
## 3 0.000999001 0.001152693
## 4 0.000999001 0.001152693
## 5 0.000999001 0.001152693
## 6 0.000999001 0.001152693
## 7 0.000999001 0.001152693
## 8 0.000999001 0.001152693
## 9 0.003996004 0.004281433
## 10 0.000999001 0.001152693
## 11 0.502497502 0.502497502
## 12 0.000999001 0.001152693
## 13 0.000999001 0.001152693
## 14 0.000999001 0.001152693
## 15 0.000999001 0.001152693#comparing all methods with each other
#evidence that everything but pitfall vs funnel traps capture significantly different reptile communities.6.5.2 Jaccard dissimilarity
all_j <- all %>% mutate(assessment.method = factor(assessment.method))
all.dist2 <- vegdist(rep_com_jacc_all, 'jaccard')
(all.adonis <- adonis2(all.dist2 ~ assessment.method, data=all_j))## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = all.dist2 ~ assessment.method, data = all_j)
## Df SumOfSqs R2 F Pr(>F)
## assessment.method 5 9.069 0.18829 4.7785 0.001 ***
## Residual 103 39.098 0.81171
## Total 108 48.168 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Permutation test: significant difference between assessment methods for reptile communities.
#Let's which method is different from which
(adonis.methods.jacc <- adonis.pair(all.dist2, all$assessment.method))## combination SumsOfSqs MeanSqs F.Model R2
## 1 camera <-> cover board 3.0365859 3.0365859 10.3763129 0.29331245
## 2 camera <-> funnel 2.6272329 2.6272329 7.0453213 0.18043958
## 3 camera <-> incidentals 2.3585467 2.3585467 6.5429368 0.20105551
## 4 camera <-> pitfall 2.7649371 2.7649371 7.6180333 0.19228701
## 5 camera <-> spotlighting 2.8379001 2.8379001 9.2684496 0.27859578
## 6 cover board <-> funnel 2.6169850 2.6169850 6.8014219 0.14849805
## 7 cover board <-> incidentals 1.3558707 1.3558707 3.5952112 0.09824267
## 8 cover board <-> pitfall 2.3940005 2.3940005 6.3568896 0.14015268
## 9 cover board <-> spotlighting 1.1861516 1.1861516 3.5283038 0.10218584
## 10 funnel <-> incidentals 1.0036015 1.0036015 2.3528111 0.05555265
## 11 funnel <-> pitfall 0.3502649 0.3502649 0.8346477 0.01782116
## 12 funnel <-> spotlighting 2.0596818 2.0596818 5.2045248 0.12046249
## 13 incidentals <-> pitfall 1.1301393 1.1301393 2.6998874 0.06322938
## 14 incidentals <-> spotlighting 0.5189479 0.5189479 1.3308740 0.03992917
## 15 pitfall <-> spotlighting 2.1228472 2.1228472 5.4802619 0.12604022
## P.value P.value.corrected
## 1 0.000999001 0.001152693
## 2 0.000999001 0.001152693
## 3 0.000999001 0.001152693
## 4 0.000999001 0.001152693
## 5 0.000999001 0.001152693
## 6 0.000999001 0.001152693
## 7 0.000999001 0.001152693
## 8 0.000999001 0.001152693
## 9 0.000999001 0.001152693
## 10 0.000999001 0.001152693
## 11 0.597402597 0.597402597
## 12 0.000999001 0.001152693
## 13 0.000999001 0.001152693
## 14 0.156843157 0.168046239
## 15 0.000999001 0.001152693#comparing all methods with each other
#evidence that everything but pitfall vs funnel traps capture significantly different reptile communities.7 Detection probability
Detection probability between methods for reptile families (maybe also terrestrial vs arboreal)
Binary linear logistic model for the likelihood to detect a reptile of the different families for the methods across sites.
7.1 Prepare data
#Number of individuals of each species captured by each method per survey plot for every survey day.
reptile_community4 <- reptile_data %>%
unite(site.plot, c(site, plot), sep = ".", remove = FALSE) %>%
select(-plot) %>%
select(site, site.plot, date, assessment.method, scientific.name) %>% #selected relevant grouping variables
group_by_at(1:4) %>%
mutate(number = n()) %>% #created a numbers column that counts the number of duplicate observations
ungroup() %>%
unique() %>% #merged the duplicate rows
pivot_wider(names_from = "scientific.name", #transformed dataset from long to wide format
values_from = "number",
values_fill = list(number=0)) %>% #replaced NAs with 0
arrange(site)
#Turn abundances for each species into a binary present/absent (1/0).
rep_com_binary <- reptile_community4 %>%
select(where(is.numeric)) %>%
decostand(method="pa", MARGIN = 1) %>%
cbind(reptile_community4[,c(2,4)]) %>%
select(site.plot, assessment.method, everything())
#Make data frame with each method for each site.plot
all.methods <- cbind.data.frame(site.plot=rep(unique(rep_com_binary$site.plot),each=6),asssessment.method=rep(unique(rep_com_binary$assessment.method),24))
#Add merged column
all.methods$site.plot.method <- paste(all.methods$site.plot,all.methods$asssessment.method)
#Find site.plot.methods that are missing from the rep_com_binary because they never documented any species at that plot
missing <- setdiff(all.methods$site.plot.method,paste(rep_com_binary$site.plot,rep_com_binary$assessment.method))
#Make data frame of those missing values with a zero filled in for each species
missing <- cbind.data.frame(all.methods[which(all.methods$site.plot.method %in% missing),1:2],as.data.frame(matrix(0, ncol=ncol(rep_com_binary)-2, nrow=length(missing))))
colnames(missing) <- colnames(rep_com_binary)
#Add missing entries
rep_com_binary <- rep_com_binary %>%
bind_rows(missing)
#Convert to long format
rep_com_binary.long <- pivot_longer(rep_com_binary, cols = -c(site.plot, assessment.method),
names_to = "scientific.name", values_to = "value")
#Get a count of days when a species was present for a given method for a given plot
days.present <- rep_com_binary.long %>%
group_by(site.plot,assessment.method,scientific.name) %>%
summarize(days.present = sum(value,na.rm=T))
#Add the survey days information (days surveyed)
survey_days <- read.csv("survey_days.csv")
days.present <- merge(days.present,survey_days,by="site.plot",all.x=T,all.y=T)
#Make a new column that is the total number of days a species was present per plot, summed across all methods and exclude plots where species was never found.
days.present <- days.present %>%
group_by(site.plot,scientific.name) %>%
mutate(count.per.site.plot = sum(days.present,na.rm=T),
assessment.method = factor(assessment.method, levels = c("pitfall", "funnel","incidentals","spotlighting","cover board","camera"))) %>%
filter(count.per.site.plot > 0) %>% #Remove plots where species were never found
select(-count.per.site.plot) %>%
mutate(det.prob = days.present/days) %>% #Probability of detecting a species via a given method at a given plot.
ungroup()This data frame now contains rows for species only for the plots where it was found but for each plot where it was found, it contains rows for each method, regardless whether that method detected it or not.
The days.present column shows the number of days when a
species was detected in that plot by that method, the days
column shows the number of days when that method was used at that plot
and the det.prob column shows the probability of detecting
a species via that method at that plot.
Include lifestyle information.
#Join with data frame that has plot and coordinate information
det.plot <- days.present %>%
left_join(coord)
#Lifestyles for species per plot
det.plot.groups <- det.plot %>%
mutate(assessment.method = factor(assessment.method, levels = c("pitfall", "funnel", "incidentals",
"spotlighting",
"cover board","camera"))) %>%
add_column(lifestyle = ifelse(.$scientific.name %in% c("Cryptoblepharus australis",
"Cryptoblepharus pannosus",
"Diplodactylus vittatus", "Gehyra dubia",
"Amalosia rhombifer", "Cryptoblepharus adamsi",
"Cryptoblepharus metallicus", "Oedura castelnaui",
"Cryptoblepharus virgatus", "Chlamydosaurus kingii",
"Varanus tristis", "Strophurus williamsi",
"Oedura coggeri", "Varanus scalaris",
"Egernia striolata", "Christinus marmoratus",
"Boiga irregularis", "Hoplocephalus bitorquatus"),
"Arboreal Species", "Ground-dwelling Species"))7.2 Bayesian Analysis - across all methods for all reptiles
7.2.1 Priors
Defining priors for Bayesian models.
Priors for the intercept.
det.plot.groups %>%
filter(days.present > 0) %>%
summarise(median(days.present))## # A tibble: 1 × 1
## `median(days.present)`
## <dbl>
## 1 2#2
det.plot.groups %>%
filter(days.present > 0) %>%
summarise(mad(days.present))## # A tibble: 1 × 1
## `mad(days.present)`
## <dbl>
## 1 1.48#1.5Priors for the slope.
sd(det.plot.groups$days.present)/apply(model.matrix(~assessment.method*lat, data = det.plot.groups), 2, sd)## (Intercept) assessment.methodfunnel
## Inf 7.2485665
## assessment.methodincidentals assessment.methodspotlighting
## 7.2485665 7.2485665
## assessment.methodcover board assessment.methodcamera
## 7.2485665 7.2485665
## lat assessment.methodfunnel:lat
## 0.4340015 0.3212041
## assessment.methodincidentals:lat assessment.methodspotlighting:lat
## 0.3212041 0.3212041
## assessment.methodcover board:lat assessment.methodcamera:lat
## 0.3212041 0.3212041#4.7 <- too high, let's go for 1.57.2.2 Fit model
7.2.2.1 Model 1
priors1 <- prior(normal(2,1.5), class = "Intercept") +
prior(normal(0,1.5), class = "b") +
prior(student_t(3,0,1.5), class = "sd")
methods.form1 <- bf(days.present ~ assessment.method*lat + offset(log(days)) + (1|site/site.plot) +
(1|scientific.name), family=poisson(link='log'))
methods.det.brm1 <- brm(methods.form1,
data = det.plot.groups,
prior = priors1,
sample_prior = "only",
refresh = 0,
chains = 3, cores = 3,
iter = 5000,
thin = 5,
seed = 1,
warmup = 1000,
backend = 'cmdstanr')## Running MCMC with 3 parallel chains...
##
## Chain 1 finished in 0.5 seconds.
## Chain 2 finished in 0.5 seconds.
## Chain 3 finished in 0.5 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 0.5 seconds.
## Total execution time: 0.6 seconds.7.2.2.1.1 Predictions
methods.det.brm1 %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)
methods.det.brm1 %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)
7.2.2.1.2 Prior vs Posterior
methods.det.brm1b <- methods.det.brm1 %>%
update(sample_prior = "yes",
refresh = 0,
seed = 1,
cores = 3,
backend = "cmdstanr")## Running MCMC with 3 parallel chains...
##
## Chain 3 finished in 130.0 seconds.
## Chain 1 finished in 169.1 seconds.
## Chain 2 finished in 172.6 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 157.2 seconds.
## Total execution time: 172.7 seconds.7.2.2.1.3 Predictions
methods.det.brm1b %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)
methods.det.brm1b %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)
7.2.2.2 Model 2
priors1 <- prior(normal(2,1.5), class = "Intercept") +
prior(normal(0,1.5), class = "b") +
prior(student_t(3,0,1.5), class = "sd")
methods.form2 <- bf(days.present ~ assessment.method*lat + offset(log(days)) + (1|site/site.plot) +
(1|scientific.name), family="negbinomial")
methods.det.brm2 <- brm(methods.form2,
data = det.plot.groups,
prior = priors1,
sample_prior = "only",
refresh = 0,
chains = 3, cores = 3,
iter = 5000,
thin = 5,
seed = 1,
warmup = 1000,
backend = 'cmdstanr')## Running MCMC with 3 parallel chains...
##
## Chain 3 finished in 3.0 seconds.
## Chain 2 finished in 4.3 seconds.
## Chain 1 finished in 4.5 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 4.0 seconds.
## Total execution time: 4.8 seconds.7.2.2.2.1 Predictions
methods.det.brm2 %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)
methods.det.brm2 %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)
7.2.2.2.2 Prior vs Posterior
methods.det.brm2b <- methods.det.brm2 %>%
update(sample_prior = "yes",
refresh = 0,
seed = 1,
cores = 3,
backend = "cmdstanr")## Running MCMC with 3 parallel chains...
##
## Chain 2 finished in 46.3 seconds.
## Chain 1 finished in 46.7 seconds.
## Chain 3 finished in 71.9 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 55.0 seconds.
## Total execution time: 72.1 seconds.7.2.2.2.3 Predictions
methods.det.brm2b %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)
methods.det.brm2b %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)
7.2.2.3 Model 3
priors1 <- prior(normal(2,1.5), class = "Intercept") +
prior(normal(0,1.5), class = "b") +
prior(student_t(3,0,1.5), class = "sd")
methods.form3 <- bf(days.present ~ assessment.method*lat + offset(log(days)) + (1|site/site.plot) +
(1|scientific.name), family="negbinomial2")
methods.det.brm3 <- brm(methods.form3,
data = det.plot.groups,
prior = priors1,
sample_prior = "only",
refresh = 0,
chains = 3, cores = 3,
iter = 5000,
thin = 5,
seed = 1,
warmup = 1000,
backend = 'cmdstanr')## Running MCMC with 3 parallel chains...
##
## Chain 1 finished in 0.5 seconds.
## Chain 2 finished in 0.5 seconds.
## Chain 3 finished in 0.5 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 0.5 seconds.
## Total execution time: 0.6 seconds.7.2.2.3.1 Predictions
methods.det.brm3 %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)
methods.det.brm3 %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)
7.2.2.3.2 Prior vs Posterior
methods.det.brm3b <- methods.det.brm3 %>%
update(sample_prior = "yes",
refresh = 0,
seed = 1,
cores = 3,
backend = "cmdstanr")## Running MCMC with 3 parallel chains...
##
## Chain 2 finished in 45.9 seconds.
## Chain 3 finished in 46.7 seconds.
## Chain 1 finished in 48.1 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 46.9 seconds.
## Total execution time: 48.1 seconds.7.2.2.3.3 Predictions
methods.det.brm3b %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)
methods.det.brm3b %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)
7.2.2.4 Compare Models
(m.1b <- methods.det.brm1b %>% loo())## Warning: Found 13 observations with a pareto_k > 0.7 in model '.'. It is
## recommended to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.##
## Computed from 2400 by 2070 log-likelihood matrix
##
## Estimate SE
## elpd_loo -3055.6 132.6
## p_loo 287.6 26.1
## looic 6111.2 265.3
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 2040 98.6% 245
## (0.5, 0.7] (ok) 17 0.8% 78
## (0.7, 1] (bad) 9 0.4% 27
## (1, Inf) (very bad) 4 0.2% 3
## See help('pareto-k-diagnostic') for details.(m.2b <- methods.det.brm2b %>% loo())## Warning: Found 7 observations with a pareto_k > 0.7 in model '.'. It is
## recommended to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.##
## Computed from 2400 by 2070 log-likelihood matrix
##
## Estimate SE
## elpd_loo -2150.7 62.1
## p_loo 61.7 5.7
## looic 4301.3 124.2
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 2060 99.5% 379
## (0.5, 0.7] (ok) 3 0.1% 245
## (0.7, 1] (bad) 7 0.3% 56
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.(m.3b <- methods.det.brm3b %>% loo())## Warning: Found 6 observations with a pareto_k > 0.7 in model '.'. It is
## recommended to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.##
## Computed from 2400 by 2070 log-likelihood matrix
##
## Estimate SE
## elpd_loo -2150.5 62.1
## p_loo 61.4 5.5
## looic 4301.0 124.1
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 2058 99.4% 431
## (0.5, 0.7] (ok) 6 0.3% 213
## (0.7, 1] (bad) 6 0.3% 63
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.loo_compare(loo(methods.det.brm1b), loo(methods.det.brm2b), loo(methods.det.brm3b))## Warning: Found 13 observations with a pareto_k > 0.7 in model
## 'methods.det.brm1b'. It is recommended to set 'moment_match = TRUE' in order to
## perform moment matching for problematic observations.## Warning: Found 7 observations with a pareto_k > 0.7 in model
## 'methods.det.brm2b'. It is recommended to set 'moment_match = TRUE' in order to
## perform moment matching for problematic observations.## Warning: Found 6 observations with a pareto_k > 0.7 in model
## 'methods.det.brm3b'. It is recommended to set 'moment_match = TRUE' in order to
## perform moment matching for problematic observations.## elpd_diff se_diff
## methods.det.brm3b 0.0 0.0
## methods.det.brm2b -0.2 0.4
## methods.det.brm1b -905.1 94.1Model methods.det.brm2b was selected as best model based
on loo estimates.
7.2.3 Diagnostics
7.2.3.1 Trace plots
methods.det.brm2b$fit %>% stan_trace()
#### Autocorrelation plots
methods.det.brm2b$fit %>% stan_ac()
#### Rhat statistic
methods.det.brm2b$fit %>% stan_rhat()
#### Effective sampling size
methods.det.brm2b$fit %>% stan_ess()
#### Posterior predictive check plot
methods.det.brm2b %>% pp_check(x="lat", type="intervals")
#### DHARMa residuals
set.seed(5)
preds <- posterior_predict(methods.det.brm2b, ndraws=250, summary=FALSE)
method.resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = det.plot.groups$days.present,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE)
method.resids %>% plot()
#### Dispersion test
method.resids %>% testDispersion()
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.38588, p-value = 0.048
## alternative hypothesis: two.sided7.2.3.2 Zero-inflation test
method.resids %>% testZeroInflation()
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.99452, p-value = 0.744
## alternative hypothesis: two.sided7.2.4 Investigation
7.2.4.1 Methods pairwise comparison for the average latitude
(diff.det.meth.avg <- methods.det.brm2b %>%
emmeans(~assessment.method|lat) %>%
regrid(transform = "none") %>%
pairs() %>%
gather_emmeans_draws() %>%
mutate(f.change = exp(.value)) %>%
summarise("Average fractional change" = median(f.change),
"Lower HDI" = HDInterval::hdi(f.change)[1],
"Upper HDI" = HDInterval::hdi(f.change)[2],
"Probability of difference" = sum(.value > 0)/n()) %>%
select(-lat))## # A tibble: 15 × 5
## # Groups: contrast [15]
## contrast Average fractional chang…¹ `Lower HDI` `Upper HDI`
## <fct> <dbl> <dbl> <dbl>
## 1 pitfall - funnel 1.11 0.832 1.47
## 2 pitfall - incidentals 6.36 4.37 8.76
## 3 pitfall - spotlighting 3.99 2.69 5.45
## 4 pitfall - cover board 5.39 3.68 7.52
## 5 pitfall - camera 27.2 15.4 42.3
## 6 funnel - incidentals 5.74 3.91 7.89
## 7 funnel - spotlighting 3.61 2.48 5.05
## 8 funnel - cover board 4.88 3.10 6.71
## 9 funnel - camera 24.4 14.2 39.3
## 10 incidentals - spotlighting 0.630 0.418 0.889
## 11 incidentals - cover board 0.848 0.544 1.19
## 12 incidentals - camera 4.27 2.36 6.78
## 13 spotlighting - cover board 1.35 0.900 1.91
## 14 spotlighting - camera 6.81 3.64 10.7
## 15 cover board - camera 5.05 2.94 8.04
## # ℹ abbreviated name: ¹`Average fractional change`
## # ℹ 1 more variable: `Probability of difference` <dbl>7.3 Bayesian Analysis - across all methods for arboreal vs ground-dwelling reptiles
7.3.1 Fit model
7.3.1.1 Model 1
priors1 <- prior(normal(2,1.5), class = "Intercept") +
prior(normal(0,1.5), class = "b") +
prior(student_t(3,0,1.5), class = "sd")
methods.form1 <- bf(days.present ~ assessment.method*lifestyle + offset(log(days)) + (1|site/site.plot) +
(1|scientific.name), family="negbinomial")
methods.det.life.brm1 <- brm(methods.form1,
data = det.plot.groups,
prior = priors1,
sample_prior = "only",
refresh = 0,
chains = 3, cores = 3,
iter = 5000,
thin = 5,
seed = 8,
warmup = 1000,
backend = 'cmdstanr')## Running MCMC with 3 parallel chains...
##
## Chain 1 finished in 3.1 seconds.
## Chain 3 finished in 3.8 seconds.
## Chain 2 finished in 3.9 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 3.6 seconds.
## Total execution time: 4.0 seconds.7.3.1.1.1 Predictions
methods.det.life.brm1 %>% ggpredict(~assessment.method|lifestyle) %>% plot(add.data = TRUE)
methods.det.life.brm1 %>% ggpredict(~lifestyle|assessment.method) %>% plot(add.data = TRUE)
7.3.1.1.2 Prior vs Posterior
methods.det.life.brm1b <- methods.det.life.brm1 %>%
update(sample_prior = "yes",
refresh = 0,
seed = 8,
cores = 3,
backend = "cmdstanr")## Running MCMC with 3 parallel chains...
##
## Chain 2 finished in 34.8 seconds.
## Chain 1 finished in 35.3 seconds.
## Chain 3 finished in 35.5 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 35.2 seconds.
## Total execution time: 35.6 seconds.7.3.1.1.3 Predictions
methods.det.life.brm1b %>% ggpredict(~assessment.method|lifestyle) %>% plot(add.data = TRUE)
methods.det.life.brm1b %>% ggpredict(~lifestyle|assessment.method) %>% plot(add.data = TRUE)
7.3.1.2 Model 2
priors2 <- prior(normal(2,1.5), class = "Intercept") +
prior(normal(0,1.5), class = "b") +
prior(student_t(3,0,1.5), class = "sd")
methods.form2 <- bf(days.present ~ assessment.method*lifestyle + offset(log(days)) + (1|site/site.plot) +
(1|scientific.name), family="negbinomial2")
methods.det.life.brm2 <- brm(methods.form2,
data = det.plot.groups,
prior = priors2,
sample_prior = "only",
refresh = 0,
chains = 3, cores = 3,
iter = 5000,
thin = 5,
seed = 8,
warmup = 1000,
backend = 'cmdstanr')## Running MCMC with 3 parallel chains...
##
## Chain 1 finished in 0.5 seconds.
## Chain 2 finished in 0.5 seconds.
## Chain 3 finished in 0.5 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 0.5 seconds.
## Total execution time: 0.6 seconds.7.3.1.2.1 Predictions
methods.det.life.brm2 %>% ggpredict(~assessment.method|lifestyle) %>% plot(add.data = TRUE)
methods.det.life.brm2 %>% ggpredict(~lifestyle|assessment.method) %>% plot(add.data = TRUE)
7.3.1.2.2 Prior vs Posterior
methods.det.life.brm2b <- methods.det.life.brm2 %>%
update(sample_prior = "yes",
refresh = 0,
seed = 8,
cores = 3,
backend = "cmdstanr")## Running MCMC with 3 parallel chains...
##
## Chain 2 finished in 34.6 seconds.
## Chain 3 finished in 36.3 seconds.
## Chain 1 finished in 59.4 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 43.4 seconds.
## Total execution time: 59.5 seconds.7.3.1.2.3 Predictions
methods.det.life.brm2b %>% ggpredict(~assessment.method|lifestyle) %>% plot(add.data = TRUE)
methods.det.life.brm2b %>% ggpredict(~lifestyle|assessment.method) %>% plot(add.data = TRUE)
7.3.1.3 Compare Models
(m.1b <- methods.det.life.brm1b %>% loo())## Warning: Found 2 observations with a pareto_k > 0.7 in model '.'. It is
## recommended to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.##
## Computed from 2400 by 2070 log-likelihood matrix
##
## Estimate SE
## elpd_loo -1993.7 58.4
## p_loo 66.9 5.8
## looic 3987.4 116.9
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 2060 99.5% 257
## (0.5, 0.7] (ok) 8 0.4% 31
## (0.7, 1] (bad) 1 0.0% 29
## (1, Inf) (very bad) 1 0.0% 13
## See help('pareto-k-diagnostic') for details.(m.2b <- methods.det.life.brm2b %>% loo())## Warning: Found 3 observations with a pareto_k > 0.7 in model '.'. It is
## recommended to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.##
## Computed from 2400 by 2070 log-likelihood matrix
##
## Estimate SE
## elpd_loo -1993.1 58.4
## p_loo 65.9 5.5
## looic 3986.3 116.7
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 2062 99.6% 355
## (0.5, 0.7] (ok) 5 0.2% 196
## (0.7, 1] (bad) 2 0.1% 65
## (1, Inf) (very bad) 1 0.0% 36
## See help('pareto-k-diagnostic') for details.loo_compare(loo(methods.det.life.brm1b), loo(methods.det.life.brm2b))## Warning: Found 2 observations with a pareto_k > 0.7 in model
## 'methods.det.life.brm1b'. It is recommended to set 'moment_match = TRUE' in
## order to perform moment matching for problematic observations.## Warning: Found 3 observations with a pareto_k > 0.7 in model
## 'methods.det.life.brm2b'. It is recommended to set 'moment_match = TRUE' in
## order to perform moment matching for problematic observations.## elpd_diff se_diff
## methods.det.life.brm2b 0.0 0.0
## methods.det.life.brm1b -0.5 0.6Model methods.det.life.brm1b was selected as best model
based on loo estimates.
7.3.2 Diagnostics
7.3.2.1 Trace plots
methods.det.life.brm1b$fit %>% stan_trace()
#### Autocorrelation plots
methods.det.life.brm1b$fit %>% stan_ac()
#### Rhat statistic
methods.det.life.brm1b$fit %>% stan_rhat()
#### Effective sampling size
methods.det.life.brm1b$fit %>% stan_ess()
7.3.2.2 DHARMa residuals
set.seed(5)
preds <- posterior_predict(methods.det.life.brm2b, ndraws=250, summary=FALSE)
method.resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = det.plot.groups$days.present,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE)
method.resids %>% plot()
#### Dispersion test
method.resids %>% testDispersion()
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.83255, p-value = 0.656
## alternative hypothesis: two.sided7.3.2.3 Zero-inflation test
method.resids %>% testZeroInflation()
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0084, p-value = 0.648
## alternative hypothesis: two.sided7.3.3 Investigation
(diff.det.life.meth.rev <- methods.det.life.brm2b %>%
emmeans(~assessment.method|lifestyle) %>%
regrid(transform = "none") %>%
pairs() %>%
#pairs(reverse = TRUE) %>% to reverse the contrasts
gather_emmeans_draws() %>%
mutate(f.change = exp(.value)) %>%
summarise("Average fractional change" = median(f.change),
"Lower HDI" = HDInterval::hdi(f.change)[1],
"Upper HDI" = HDInterval::hdi(f.change)[2],
"Probability of difference" = sum(.value > 0)/n()) %>%
arrange(lifestyle))## # A tibble: 30 × 6
## # Groups: contrast [15]
## contrast lifestyle Average fractional c…¹ `Lower HDI` `Upper HDI`
## <fct> <fct> <dbl> <dbl> <dbl>
## 1 pitfall - funnel Arboreal… 3.85 1.70 6.62
## 2 pitfall - incidenta… Arboreal… 2.98 1.52 4.95
## 3 pitfall - spotlight… Arboreal… 0.541 0.309 0.864
## 4 pitfall - cover boa… Arboreal… 0.452 0.279 0.718
## 5 pitfall - camera Arboreal… 27.9 7.01 71.6
## 6 funnel - incidentals Arboreal… 0.770 0.335 1.44
## 7 funnel - spotlighti… Arboreal… 0.139 0.0643 0.241
## 8 funnel - cover board Arboreal… 0.118 0.0551 0.199
## 9 funnel - camera Arboreal… 7.26 1.54 20.6
## 10 incidentals - spotl… Arboreal… 0.180 0.0984 0.293
## # ℹ 20 more rows
## # ℹ abbreviated name: ¹`Average fractional change`
## # ℹ 1 more variable: `Probability of difference` <dbl>7.4 Visualisation
pal <- c("#FF8300", "#D14103","#0CC170","black","#4E84C4","#8348D8")
(det.life.plot <- ggplot(det.plot.groups, aes(x=assessment.method,y=det.prob)) +
geom_boxplot(aes(color = assessment.method,
color = after_scale(darken(color, .1, space = "HLS")),
fill = after_scale(desaturate(lighten(color, .8), .4))),
width = .6, outlier.shape = NA) +
geom_point(aes(color = assessment.method, color = after_scale(darken(color, .1, space = "HLS"))), fill = "white", shape = 21, stroke = .4, size = 3.5, side = "1",
transformation = PositionIdentity) +
geom_point(aes(fill = assessment.method), color = "transparent", shape = 21, stroke = .4, size = 3.5, alpha = .3, side = "1", transformation = PositionIdentity) +
facet_wrap(~lifestyle) +
scale_y_continuous(labels = scales::percent,
breaks = seq(0, 0.9, by = 0.1),
limits = c(0,0.9),
name = "Detection Probability") +
scale_x_discrete(name = "", guide = "none") +
scale_fill_manual(values = pal, guide = "none") +
scale_colour_manual(values = pal, name = "", labels = c("Pitfall Trap","Funnel Trap","Incidental","Spotlighting","Arboreal Cover Board","Camera Trap")) +
my.theme() +
theme(panel.grid.major.y = element_line(colour = "grey", linewidth = 0.2, linetype = "dotted"),
panel.grid.minor.y = element_line(colour = "grey", linewidth = 0.1, linetype = "dotted"),
legend.position = c(0.5, -0.07),
legend.background = element_rect(fill = "transparent"),
legend.direction = "horizontal",
plot.margin = unit(c(5, 10, 12, 10), units = "mm"),
panel.border = element_rect(fill = NA, color = "black"),
strip.background = element_rect(fill = "lightgrey")))## Warning: Duplicated aesthetics after name standardisation: colour
## Duplicated aesthetics after name standardisation: colour## Warning in geom_point(aes(color = assessment.method, color =
## after_scale(darken(color, : Ignoring unknown parameters: `side` and
## `transformation`## Warning in geom_point(aes(fill = assessment.method), color = "transparent", :
## Ignoring unknown parameters: `side` and `transformation`
8 Session Info
sessionInfo()## R version 4.3.0 (2023-04-21)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Ventura 13.3.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Australia/Brisbane
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] report_0.5.7 HDInterval_0.2.4 patchwork_1.1.2
## [4] EcolUtils_0.1 vegan_2.6-4 lattice_0.21-8
## [7] permute_0.9-7 ggridges_0.5.4 tidybayes_3.0.4
## [10] emmeans_1.8.6 DHARMa_0.4.6 ggeffects_1.2.2
## [13] rstan_2.21.8 StanHeaders_2.21.0-7 brms_2.19.0
## [16] Rcpp_1.0.10 cmdstanr_0.5.3 colorspace_2.1-0
## [19] lubridate_1.9.2 forcats_1.0.0 stringr_1.5.0
## [22] dplyr_1.1.2 purrr_1.0.1 readr_2.1.4
## [25] tidyr_1.3.0 tibble_3.2.1 ggplot2_3.4.2
## [28] tidyverse_2.0.0 knitr_1.42
##
## loaded via a namespace (and not attached):
## [1] tensorA_0.36.2 rstudioapi_0.14 jsonlite_1.8.4
## [4] magrittr_2.0.3 estimability_1.4.1 farver_2.1.1
## [7] nloptr_2.0.3 rmarkdown_2.21 vctrs_0.6.2
## [10] minqa_1.2.5 base64enc_0.1-3 htmltools_0.5.5
## [13] haven_2.5.2 distributional_0.3.2 sass_0.4.6
## [16] bslib_0.4.2 htmlwidgets_1.6.2 plyr_1.8.8
## [19] zoo_1.8-12 cachem_1.0.8 igraph_1.4.2
## [22] iterators_1.0.14 mime_0.12 lifecycle_1.0.3
## [25] pkgconfig_2.0.3 gap_1.5-1 colourpicker_1.2.0
## [28] Matrix_1.5-4 R6_2.5.1 fastmap_1.1.1
## [31] shiny_1.7.4 digest_0.6.31 ps_1.7.5
## [34] qgam_1.3.4 crosstalk_1.2.0 labeling_0.4.2
## [37] fansi_1.0.4 timechange_0.2.0 abind_1.4-5
## [40] mgcv_1.8-42 compiler_4.3.0 doParallel_1.0.17
## [43] withr_2.5.0 backports_1.4.1 inline_0.3.19
## [46] shinystan_2.6.0 highr_0.10 pkgbuild_1.4.0
## [49] MASS_7.3-58.4 gtools_3.9.4 loo_2.6.0
## [52] tools_4.3.0 httpuv_1.6.11 threejs_0.3.3
## [55] glue_1.6.2 callr_3.7.3 nlme_3.1-162
## [58] promises_1.2.0.1 grid_4.3.0 checkmate_2.2.0
## [61] cluster_2.1.4 reshape2_1.4.4 generics_0.1.3
## [64] gtable_0.3.3 tzdb_0.4.0 data.table_1.14.8
## [67] hms_1.1.3 utf8_1.2.3 foreach_1.5.2
## [70] pillar_1.9.0 ggdist_3.3.0 markdown_1.7
## [73] posterior_1.4.1 later_1.3.1 splines_4.3.0
## [76] tidyselect_1.2.0 miniUI_0.1.1.1 arrayhelpers_1.1-0
## [79] gridExtra_2.3 stats4_4.3.0 xfun_0.39
## [82] bridgesampling_1.1-2 matrixStats_0.63.0 DT_0.27
## [85] stringi_1.7.12 yaml_2.3.7 boot_1.3-28.1
## [88] evaluate_0.21 codetools_0.2-19 cli_3.6.1
## [91] RcppParallel_5.1.7 shinythemes_1.2.0 xtable_1.8-4
## [94] munsell_0.5.0 processx_3.8.1 jquerylib_0.1.4
## [97] coda_0.19-4 svUnit_1.0.6 parallel_4.3.0
## [100] rstantools_2.3.1 ellipsis_0.3.2 prettyunits_1.1.1
## [103] gap.datasets_0.0.5 dygraphs_1.1.1.6 bayesplot_1.10.0
## [106] Brobdingnag_1.2-9 lme4_1.1-33 mvtnorm_1.1-3
## [109] scales_1.2.1 xts_0.13.1 insight_0.19.1
## [112] crayon_1.5.2 rlang_1.1.1 shinyjs_2.1.09 Cite Packages
cite_packages()## - Bürkner P (2017). "brms: An R Package for Bayesian Multilevel Models Using Stan." _Journal of Statistical Software_, *80*(1), 1-28. doi:10.18637/jss.v080.i01 <https://doi.org/10.18637/jss.v080.i01>. Bürkner P (2018). "Advanced Bayesian Multilevel Modeling with the R Package brms." _The R Journal_, *10*(1), 395-411. doi:10.32614/RJ-2018-017 <https://doi.org/10.32614/RJ-2018-017>. Bürkner P (2021). "Bayesian Item Response Modeling in R with brms and Stan." _Journal of Statistical Software_, *100*(5), 1-54. doi:10.18637/jss.v100.i05 <https://doi.org/10.18637/jss.v100.i05>.
## - Eddelbuettel D, François R (2011). "Rcpp: Seamless R and C++ Integration." _Journal of Statistical Software_, *40*(8), 1-18. doi:10.18637/jss.v040.i08 <https://doi.org/10.18637/jss.v040.i08>. Eddelbuettel D (2013). _Seamless R and C++ Integration with Rcpp_. Springer, New York. doi:10.1007/978-1-4614-6868-4 <https://doi.org/10.1007/978-1-4614-6868-4>, ISBN 978-1-4614-6867-7. Eddelbuettel D, Balamuta JJ (2018). "Extending extitR with extitC++: A Brief Introduction to extitRcpp." _The American Statistician_, *72*(1), 28-36. doi:10.1080/00031305.2017.1375990 <https://doi.org/10.1080/00031305.2017.1375990>.
## - Gabry J, Češnovar R (2022). _cmdstanr: R Interface to 'CmdStan'_. https://mc-stan.org/cmdstanr/, https://discourse.mc-stan.org.
## - Grolemund G, Wickham H (2011). "Dates and Times Made Easy with lubridate." _Journal of Statistical Software_, *40*(3), 1-25. <https://www.jstatsoft.org/v40/i03/>.
## - Hartig F (2022). _DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models_. R package version 0.4.6, <https://CRAN.R-project.org/package=DHARMa>.
## - Kay M (2023). _tidybayes: Tidy Data and Geoms for Bayesian Models_. doi:10.5281/zenodo.1308151 <https://doi.org/10.5281/zenodo.1308151>, R package version 3.0.4, <http://mjskay.github.io/tidybayes/>.
## - Lenth R (2023). _emmeans: Estimated Marginal Means, aka Least-Squares Means_. R package version 1.8.6, <https://CRAN.R-project.org/package=emmeans>.
## - Lüdecke D (2018). "ggeffects: Tidy Data Frames of Marginal Effects from Regression Models." _Journal of Open Source Software_, *3*(26), 772. doi:10.21105/joss.00772 <https://doi.org/10.21105/joss.00772>.
## - Makowski D, Lüdecke D, Patil I, Thériault R, Ben-Shachar M, Wiernik B (2023). "Automated Results Reporting as a Practical Tool to Improve Reproducibility and Methodological Best Practices Adoption." _CRAN_. <https://easystats.github.io/report/>.
## - Meredith M, Kruschke J (2022). _HDInterval: Highest (Posterior) Density Intervals_. R package version 0.2.4, <https://CRAN.R-project.org/package=HDInterval>.
## - Müller K, Wickham H (2023). _tibble: Simple Data Frames_. R package version 3.2.1, <https://CRAN.R-project.org/package=tibble>.
## - Oksanen J, Simpson G, Blanchet F, Kindt R, Legendre P, Minchin P, O'Hara R, Solymos P, Stevens M, Szoecs E, Wagner H, Barbour M, Bedward M, Bolker B, Borcard D, Carvalho G, Chirico M, De Caceres M, Durand S, Evangelista H, FitzJohn R, Friendly M, Furneaux B, Hannigan G, Hill M, Lahti L, McGlinn D, Ouellette M, Ribeiro Cunha E, Smith T, Stier A, Ter Braak C, Weedon J (2022). _vegan: Community Ecology Package_. R package version 2.6-4, <https://CRAN.R-project.org/package=vegan>.
## - Pedersen T (2022). _patchwork: The Composer of Plots_. R package version 1.1.2, <https://CRAN.R-project.org/package=patchwork>.
## - R Core Team (2023). _R: A Language and Environment for Statistical Computing_. R Foundation for Statistical Computing, Vienna, Austria. <https://www.R-project.org/>.
## - Salazar G (2023). _EcolUtils: Utilities for community ecology analysis_. R package version 0.1, <https://github.com/GuillemSalazar/EcolUtils>.
## - Sarkar D (2008). _Lattice: Multivariate Data Visualization with R_. Springer, New York. ISBN 978-0-387-75968-5, <http://lmdvr.r-forge.r-project.org>.
## - Simpson G (2022). _permute: Functions for Generating Restricted Permutations of Data_. R package version 0.9-7, <https://CRAN.R-project.org/package=permute>.
## - Stan Development Team (2020). "StanHeaders: Headers for the R interface to Stan." R package version 2.21.0-6, <https://mc-stan.org/>.
## - Stan Development Team (2023). "RStan: the R interface to Stan." R package version 2.21.8, <https://mc-stan.org/>.
## - Wickham H (2016). _ggplot2: Elegant Graphics for Data Analysis_. Springer-Verlag New York. ISBN 978-3-319-24277-4, <https://ggplot2.tidyverse.org>.
## - Wickham H (2022). _stringr: Simple, Consistent Wrappers for Common String Operations_. R package version 1.5.0, <https://CRAN.R-project.org/package=stringr>.
## - Wickham H (2023). _forcats: Tools for Working with Categorical Variables (Factors)_. R package version 1.0.0, <https://CRAN.R-project.org/package=forcats>.
## - Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). "Welcome to the tidyverse." _Journal of Open Source Software_, *4*(43), 1686. doi:10.21105/joss.01686 <https://doi.org/10.21105/joss.01686>.
## - Wickham H, François R, Henry L, Müller K, Vaughan D (2023). _dplyr: A Grammar of Data Manipulation_. R package version 1.1.2, <https://CRAN.R-project.org/package=dplyr>.
## - Wickham H, Henry L (2023). _purrr: Functional Programming Tools_. R package version 1.0.1, <https://CRAN.R-project.org/package=purrr>.
## - Wickham H, Hester J, Bryan J (2023). _readr: Read Rectangular Text Data_. R package version 2.1.4, <https://CRAN.R-project.org/package=readr>.
## - Wickham H, Vaughan D, Girlich M (2023). _tidyr: Tidy Messy Data_. R package version 1.3.0, <https://CRAN.R-project.org/package=tidyr>.
## - Wilke C (2022). _ggridges: Ridgeline Plots in 'ggplot2'_. R package version 0.5.4, <https://CRAN.R-project.org/package=ggridges>.
## - Xie Y (2023). _knitr: A General-Purpose Package for Dynamic Report Generation in R_. R package version 1.42, <https://yihui.org/knitr/>. Xie Y (2015). _Dynamic Documents with R and knitr_, 2nd edition. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 978-1498716963, <https://yihui.org/knitr/>. Xie Y (2014). "knitr: A Comprehensive Tool for Reproducible Research in R." In Stodden V, Leisch F, Peng RD (eds.), _Implementing Reproducible Computational Research_. Chapman and Hall/CRC. ISBN 978-1466561595.
## - Zeileis A, Fisher JC, Hornik K, Ihaka R, McWhite CD, Murrell P, Stauffer R, Wilke CO (2020). "colorspace: A Toolbox for Manipulating and Assessing Colors and Palettes." _Journal of Statistical Software_, *96*(1), 1-49. doi:10.18637/jss.v096.i01 <https://doi.org/10.18637/jss.v096.i01>. Zeileis A, Hornik K, Murrell P (2009). "Escaping RGBland: Selecting Colors for Statistical Graphics." _Computational Statistics & Data Analysis_, *53*(9), 3259-3270. doi:10.1016/j.csda.2008.11.033 <https://doi.org/10.1016/j.csda.2008.11.033>. Stauffer R, Mayr GJ, Dabernig M, Zeileis A (2009). "Somewhere over the Rainbow: How to Make Effective Use of Colors in Meteorological Visualizations." _Bulletin of the American Meteorological Society_, *96*(2), 203-216. doi:10.1175/BAMS-D-13-00155.1 <https://doi.org/10.1175/BAMS-D-13-00155.1>.